home *** CD-ROM | disk | FTP | other *** search
/ IRIX 6.2 Development Libraries / SGI IRIX 6.2 Development Libraries.iso / dist / complib.idb / usr / share / catman / p_man / cat3 / complib / dlatbs.z / dlatbs
Text File  |  1996-03-14  |  8KB  |  265 lines
  1.  
  2.  
  3.  
  4. DDDDLLLLAAAATTTTBBBBSSSS((((3333FFFF))))                                                          DDDDLLLLAAAATTTTBBBBSSSS((((3333FFFF))))
  5.  
  6.  
  7.  
  8. NNNNAAAAMMMMEEEE
  9.      DLATBS - solve one of the triangular systems   A *x = s*b or A'*x = s*b
  10.      with scaling to prevent overflow, where A is an upper or lower triangular
  11.      band matrix
  12.  
  13. SSSSYYYYNNNNOOOOPPPPSSSSIIIISSSS
  14.      SUBROUTINE DLATBS( UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X, SCALE,
  15.                         CNORM, INFO )
  16.  
  17.          CHARACTER      DIAG, NORMIN, TRANS, UPLO
  18.  
  19.          INTEGER        INFO, KD, LDAB, N
  20.  
  21.          DOUBLE         PRECISION SCALE
  22.  
  23.          DOUBLE         PRECISION AB( LDAB, * ), CNORM( * ), X( * )
  24.  
  25. PPPPUUUURRRRPPPPOOOOSSSSEEEE
  26.      DLATBS solves one of the triangular systems are n-element vectors, and s
  27.      is a scaling factor, usually less than or equal to 1, chosen so that the
  28.      components of x will be less than the overflow threshold.  If the
  29.      unscaled problem will not cause overflow, the Level 2 BLAS routine DTBSV
  30.      is called.  If the matrix A is singular (A(j,j) = 0 for some j), then s
  31.      is set to 0 and a non-trivial solution to A*x = 0 is returned.
  32.  
  33.  
  34. AAAARRRRGGGGUUUUMMMMEEEENNNNTTTTSSSS
  35.      UPLO    (input) CHARACTER*1
  36.              Specifies whether the matrix A is upper or lower triangular.  =
  37.              'U':  Upper triangular
  38.              = 'L':  Lower triangular
  39.  
  40.      TRANS   (input) CHARACTER*1
  41.              Specifies the operation applied to A.  = 'N':  Solve A * x = s*b
  42.              (No transpose)
  43.              = 'T':  Solve A'* x = s*b  (Transpose)
  44.              = 'C':  Solve A'* x = s*b  (Conjugate transpose = Transpose)
  45.  
  46.      DIAG    (input) CHARACTER*1
  47.              Specifies whether or not the matrix A is unit triangular.  = 'N':
  48.              Non-unit triangular
  49.              = 'U':  Unit triangular
  50.  
  51.      NORMIN  (input) CHARACTER*1
  52.              Specifies whether CNORM has been set or not.  = 'Y':  CNORM
  53.              contains the column norms on entry
  54.              = 'N':  CNORM is not set on entry.  On exit, the norms will be
  55.              computed and stored in CNORM.
  56.  
  57.      N       (input) INTEGER
  58.              The order of the matrix A.  N >= 0.
  59.  
  60.  
  61.  
  62.  
  63.                                                                         PPPPaaaaggggeeee 1111
  64.  
  65.  
  66.  
  67.  
  68.  
  69.  
  70. DDDDLLLLAAAATTTTBBBBSSSS((((3333FFFF))))                                                          DDDDLLLLAAAATTTTBBBBSSSS((((3333FFFF))))
  71.  
  72.  
  73.  
  74.      KD      (input) INTEGER
  75.              The number of subdiagonals or superdiagonals in the triangular
  76.              matrix A.  KD >= 0.
  77.  
  78.      AB      (input) DOUBLE PRECISION array, dimension (LDAB,N)
  79.              The upper or lower triangular band matrix A, stored in the first
  80.              KD+1 rows of the array. The j-th column of A is stored in the j-
  81.              th column of the array AB as follows:  if UPLO = 'U', AB(kd+1+i-
  82.              j,j) = A(i,j) for max(1,j-kd)<=i<=j; if UPLO = 'L', AB(1+i-j,j)
  83.              = A(i,j) for j<=i<=min(n,j+kd).
  84.  
  85.      LDAB    (input) INTEGER
  86.              The leading dimension of the array AB.  LDAB >= KD+1.
  87.  
  88.      X       (input/output) DOUBLE PRECISION array, dimension (N)
  89.              On entry, the right hand side b of the triangular system.  On
  90.              exit, X is overwritten by the solution vector x.
  91.  
  92.      SCALE   (output) DOUBLE PRECISION
  93.              The scaling factor s for the triangular system A * x = s*b  or
  94.              A'* x = s*b.  If SCALE = 0, the matrix A is singular or badly
  95.              scaled, and the vector x is an exact or approximate solution to
  96.              A*x = 0.
  97.  
  98.      CNORM   (input or output) DOUBLE PRECISION array, dimension (N)
  99.  
  100.              If NORMIN = 'Y', CNORM is an input argument and CNORM(j) contains
  101.              the norm of the off-diagonal part of the j-th column of A.  If
  102.              TRANS = 'N', CNORM(j) must be greater than or equal to the
  103.              infinity-norm, and if TRANS = 'T' or 'C', CNORM(j) must be
  104.              greater than or equal to the 1-norm.
  105.  
  106.              If NORMIN = 'N', CNORM is an output argument and CNORM(j) returns
  107.              the 1-norm of the offdiagonal part of the j-th column of A.
  108.  
  109.      INFO    (output) INTEGER
  110.              = 0:  successful exit
  111.              < 0:  if INFO = -k, the k-th argument had an illegal value
  112.  
  113. FFFFUUUURRRRTTTTHHHHEEEERRRR DDDDEEEETTTTAAAAIIIILLLLSSSS
  114.      A rough bound on x is computed; if that is less than overflow, DTBSV is
  115.      called, otherwise, specific code is used which checks for possible
  116.      overflow or divide-by-zero at every operation.
  117.  
  118.      A columnwise scheme is used for solving A*x = b.  The basic algorithm if
  119.      A is lower triangular is
  120.  
  121.           x[1:n] := b[1:n]
  122.           for j = 1, ..., n
  123.                x(j) := x(j) / A(j,j)
  124.                x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
  125.           end
  126.  
  127.  
  128.  
  129.                                                                         PPPPaaaaggggeeee 2222
  130.  
  131.  
  132.  
  133.  
  134.  
  135.  
  136. DDDDLLLLAAAATTTTBBBBSSSS((((3333FFFF))))                                                          DDDDLLLLAAAATTTTBBBBSSSS((((3333FFFF))))
  137.  
  138.  
  139.  
  140.      Define bounds on the components of x after j iterations of the loop:
  141.         M(j) = bound on x[1:j]
  142.         G(j) = bound on x[j+1:n]
  143.      Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.
  144.  
  145.      Then for iteration j+1 we have
  146.         M(j+1) <= G(j) / | A(j+1,j+1) |
  147.         G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
  148.                <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )
  149.  
  150.      where CNORM(j+1) is greater than or equal to the infinity-norm of column
  151.      j+1 of A, not counting the diagonal.  Hence
  152.  
  153.         G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
  154.                      1<=i<=j
  155.      and
  156.  
  157.         |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
  158.                                       1<=i< j
  159.  
  160.      Since |x(j)| <= M(j), we use the Level 2 BLAS routine DTBSV if the
  161.      reciprocal of the largest M(j), j=1,..,n, is larger than
  162.      max(underflow, 1/overflow).
  163.  
  164.      The bound on x(j) is also used to determine when a step in the columnwise
  165.      method can be performed without fear of overflow.  If the computed bound
  166.      is greater than a large constant, x is scaled to prevent overflow, but if
  167.      the bound overflows, x is set to 0, x(j) to 1, and scale to 0, and a
  168.      non-trivial solution to A*x = 0 is found.
  169.  
  170.      Similarly, a row-wise scheme is used to solve A'*x = b.  The basic
  171.      algorithm for A upper triangular is
  172.  
  173.           for j = 1, ..., n
  174.                x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j)
  175.           end
  176.  
  177.      We simultaneously compute two bounds
  178.           G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j
  179.           M(j) = bound on x(i), 1<=i<=j
  180.  
  181.      The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we add
  182.      the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.  Then the
  183.      bound on x(j) is
  184.  
  185.           M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |
  186.  
  187.                <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
  188.                          1<=i<=j
  189.  
  190.      and we can safely call DTBSV if 1/M(n) and 1/G(n) are both greater than
  191.      max(underflow, 1/overflow).
  192.  
  193.  
  194.  
  195.                                                                         PPPPaaaaggggeeee 3333
  196.  
  197.  
  198.  
  199.  
  200.  
  201.  
  202. DDDDLLLLAAAATTTTBBBBSSSS((((3333FFFF))))                                                          DDDDLLLLAAAATTTTBBBBSSSS((((3333FFFF))))
  203.  
  204.  
  205.  
  206.  
  207.  
  208.  
  209.  
  210.  
  211.  
  212.  
  213.  
  214.  
  215.  
  216.  
  217.  
  218.  
  219.  
  220.  
  221.  
  222.  
  223.  
  224.  
  225.  
  226.  
  227.  
  228.  
  229.  
  230.  
  231.  
  232.  
  233.  
  234.  
  235.  
  236.  
  237.  
  238.  
  239.  
  240.  
  241.  
  242.  
  243.  
  244.  
  245.  
  246.  
  247.  
  248.  
  249.  
  250.  
  251.  
  252.  
  253.  
  254.  
  255.  
  256.  
  257.  
  258.                                                                         PPPPaaaaggggeeee 4444
  259.  
  260.  
  261.  
  262.  
  263.  
  264.  
  265.